%=====================
%this file solves the T3 treatment (without TDA, new tax function)
%the new tax rate on interest income is solved such that the expected
%government revenue is the same in T3 as in T1 (TDA treatment)
%=====================
clearvars;
global  climit z JR 
clc;
close all


%=====================
%define environment values
%=====================
JR=16; %retirement period 
r=0.10; %per period pre-tax interest rate
b=1; %first we start with no discount factor
%endowment for workers
e=[67
74
79
85
89
93
96
97
98
98
97
95
93
90
86
];

eJr=51; %endowment for retirees

climit=47.4;
JN=100;
s=9/10; %1/s is the expected lenghth for a constant survival rate s

%z is the new tax rate on interest income
zhigh=0.1;
zlow=0;
critz=100;
amax=2500;

%====================
%Define state space and choice/value function
%=====================
NA=100; 
Avector=linspace(0,30,NA)';
Avector=Avector.^2; %more grids for small values
%for workers
vfunc=cell(JR,1); %the cell argument represents period
cfunc=cell(JR,1);

while critz>0.01
    z=(zhigh+zlow)/2

%=====================
%Solve the model
%=====================

%start with retirees
critv=1.0e4;
diffv=1.0e-7;
crit=0.0004;
NCini=5; %initially check five points
for j=1:JR
    vfunc{j}=zeros(NA,1);
    cfunc{j}=zeros(NA,1);
end %last period equals that of retirees
%for retirees
vrfunc=zeros(NA,1);
crfunc=zeros(NA,1);

while (critv>diffv)
    vrfunc_old=vrfunc;
          for ia=1:NA
            if ia==1
                cmin=climit;
            else
            cmin=crfunc(ia-1);
            end
            cmax=(eJr+(1+r)*Avector(ia)-tax_T3(eJr,r*Avector(ia)));
            cg=linspace(cmin,cmax,NCini);
            afuture=eJr+(1+r)*Avector(ia)-tax_T3(eJr,r*Avector(ia))-cg;
            vtemp=uc(cg)+b*s*interp1(Avector,vrfunc,min(amax,afuture),'spline');
            [vtemp_v, ctemp_g]=max(vtemp);
            cmax_v=cg(ctemp_g);
            vmax_v=vtemp_v;
            
            diff=(cmax-cmin)/4;
            while diff>crit

            diff=diff/2;
          
            if ((cmax_v>cmin))
                cg=cmax_v-diff;
                afuture=eJr+(1+r)*Avector(ia)-tax_T3(eJr,r*Avector(ia))-cg;
                vtemp=uc(cg)+b*s*interp1(Avector,vrfunc,min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                    continue
                end %if
            end %while internal
            if ((cmax_v<cmax))
                cg=cmax_v+diff;
                afuture=eJr+(1+r)*Avector(ia)-tax_T3(eJr,r*Avector(ia))-cg;
                vtemp=uc(cg)+b*s*interp1(Avector,vrfunc,min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                end %if
            end %while internal
            end %while
            vrfunc(ia)=vmax_v;
            crfunc(ia)=cmax_v;
     
           end %ia
            
critv=max(abs(vrfunc_old-vrfunc));
end
%for workers
vfunc_temp=zeros(NA,1);
cfunc_temp=zeros(NA,1);
vfunc{JR}=vrfunc;
cfunc{JR}=crfunc;
for j=JR-1:-1:1 
    vfuture=vfunc{j+1};
        for ia=1:NA
            if ia==1
                cmin=climit;
            else
            cmin=cfunc_temp(ia-1);
            end
            cmax=(e(j)+(1+r)*Avector(ia)-tax_T3(e(j),r*Avector(ia)));
            cg=linspace(cmin,cmax,NCini);
            afuture=e(j)+(1+r)*Avector(ia)-tax_T3(e(j),r*Avector(ia))-cg;
            vtemp=uc(cg)+b*interp1(Avector,vfuture,min(amax,afuture),'spline');
            [vtemp_v, ctemp_g]=max(vtemp);
            cmax_v=cg(ctemp_g);
            vmax_v=vtemp_v;
            
            diff=(cmax-cmin)/4;
            while diff>crit

            diff=diff/2;
          
            if ((cmax_v>cmin))
                cg=cmax_v-diff;
                afuture=e(j)+(1+r)*Avector(ia)-tax_T3(e(j),r*Avector(ia))-cg;
                vtemp=uc(cg)+b*interp1(Avector,vfuture,min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                    continue
                end %if
            end %while internal
            if ((cmax_v<cmax))
                cg=cmax_v+diff;
                afuture=e(j)+(1+r)*Avector(ia)-tax_T3(e(j),r*Avector(ia))-cg;
                vtemp=uc(cg)+b*interp1(Avector,vfuture,min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                end %if
            end %while internal
            end %while
            vfunc_temp(ia)=vmax_v;
            cfunc_temp(ia)=cmax_v;
     
        end %ia
        vfunc{j}=vfunc_temp;
        cfunc{j}=cfunc_temp;
end %j


%=====================
%simulate optimal decisions
%=====================
JN=100; %maximum periods
vopt=zeros(JN,1); %optimal value
copt=zeros(JN,1); %consumption
aopt=zeros(JN+1,1); %asset
taxopt=zeros(JN,1); %tax liability
taxableopt=zeros(JN,1); %tax payment
for j=1:JR-1
            vfuture=vfunc{j+1};
            cmin=climit;
            cmax=(e(j)+(1+r)*aopt(j)-tax_T3(e(j),r*aopt(j)));
             %start with three points
            cg=linspace(cmin,cmax,NCini);
            afuture=e(j)+(1+r)*aopt(j)-tax_T3(e(j),r*aopt(j))-cg;
            vtemp=uc(cg)+b*interp1(Avector,vfuture, min(amax,afuture),'spline');
            [vtemp_v, ctemp_g]=max(vtemp);
            cmax_v=cg(ctemp_g);
            vmax_v=vtemp_v;
            
            diff=(cmax-cmin)/4;
            while diff>crit

            diff=diff/2;
          
            if ((cmax_v>cmin))
                cg=cmax_v-diff;
                afuture=e(j)+(1+r)*aopt(j)-tax_T3(e(j),r*aopt(j))-cg;
                vtemp=uc(cg)+b*interp1(Avector,vfuture, min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                    continue
                end %if
            end %while internal
            if ((cmax_v<cmax))
                cg=cmax_v+diff;
                afuture=e(j)+(1+r)*aopt(j)-tax_T3(e(j),r*aopt(j))-cg;
                vtemp=uc(cg)+b*interp1(Avector,vfuture, min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                end %if
            end %while internal
            end %while
            copt(j)=min(round(cmax_v,1),floor((e(j)+(1+r)*aopt(j)-tax_T3(e(j),r*aopt(j)))*1000)/1000);          

    aopt(j+1)=e(j)+(1+r)*aopt(j)-tax_T3(e(j),r*aopt(j))-copt(j);   
    taxopt(j)=tax_T3(e(j),r*aopt(j));
    taxableopt(j)=e(j)+r*aopt(j);
end

for j=JR:JN 
            vfuture=vrfunc;
            cmin=climit;
            cmax=(eJr+(1+r)*aopt(j)-tax_T3(eJr,r*aopt(j)));
             %start with three points
            cg=linspace(cmin,cmax,NCini);
            afuture=eJr+(1+r)*aopt(j)-tax_T3(eJr,r*aopt(j))-cg;
            vtemp=uc(cg)+b*s*interp1(Avector,vfuture, min(amax,afuture),'spline');
            [vtemp_v, ctemp_g]=max(vtemp);
            cmax_v=cg(ctemp_g);
            vmax_v=vtemp_v;
            
            diff=(cmax-cmin)/4;
            while diff>crit

            diff=diff/2;
          
            if ((cmax_v>cmin))
                cg=cmax_v-diff;
                afuture=eJr+(1+r)*aopt(j)-tax_T3(eJr,r*aopt(j))-cg;
                vtemp=uc(cg)+b*s*interp1(Avector,vfuture, min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                    continue
                end %if
            end %while internal
            if ((cmax_v<cmax))
                cg=cmax_v+diff;
                afuture=eJr+(1+r)*aopt(j)-tax_T3(eJr,r*aopt(j))-cg;
                vtemp=uc(cg)+b*s*interp1(Avector,vfuture, min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                end %if
            end %while internal
            end %while
            copt(j)=min(round(cmax_v,1),floor((eJr+(1+r)*aopt(j)-tax_T3(eJr,r*aopt(j)))*1000)/1000);          

    aopt(j+1)=eJr+(1+r)*aopt(j)-tax_T3(eJr,r*aopt(j))-copt(j);   
    taxopt(j)=tax_T3(eJr,r*aopt(j));
    taxableopt(j)=eJr+r*aopt(j);
end


%cumulative payoff if adopting the optimal strategy
paymentopt=zeros(JN,1);
j=1;
paymentopt(j)=b^(j-1)*uc(copt(j));
for j=2:JN
paymentopt(j)=paymentopt(j-1)+b^(j-1)*uc(copt(j));
end




%hand to mouth consumers
paymenthandtomouth=zeros(JN,1);
j=1;
paymenthandtomouth(j)=b^(j-1)*uc(e(j)-tax_T3(e(j),0));
for j=2:JR-1
paymenthandtomouth(j)=paymenthandtomouth(j-1)+b^(j-1)*uc(e(j)-tax_T3(e(j),0));
end
for j=JR:JN
paymenthandtomouth(j)=paymenthandtomouth(j-1)+b^(j-1)*uc(eJr-tax_T3(eJr,0));
end

chandtomouth=zeros(JN,1);
j=1;
chandtomouth(j)=e(j)-tax_T3(e(j),0);
for j=2:JR-1
chandtomouth(j)=(e(j)-tax_T3(e(j),0));
end
for j=JR:JN
chandtomouth(j)=(eJr-tax_T3(eJr,0));
end

prob_d=zeros(JN+1,1); %the probability of decease in period j
for j=JR+1:JN
    prob_d(j)=prob_d(j-1)+(1-s)*(1-prob_d(j-1));
end
prob_d(JN+1)=1; %maximum age of JN

fprintf('Expected tax revenue')
sum((1-prob_d(1:JN)).*taxopt)
critz=abs(mean(sum((1-prob_d(1:JN)).*taxopt)-211.339424739925))
%update the bound
if sum((1-prob_d(1:JN)).*taxopt)>211.339424739925
    zhigh=z;
 
else
    zlow=z;


end
end


%deal with the rounding of z in experiments
z=round(z,4);
fprintf('final tax rate')
z

for j=1:JR
    vfunc{j}=zeros(NA,1);
    cfunc{j}=zeros(NA,1);
end %last period equals that of retirees
%for retirees
vrfunc=zeros(NA,1);
crfunc=zeros(NA,1);
%=====================
%Solve the model under the final rounded tax rate z
%=====================
tic
%start with retirees
critv=1.0e4;
diffv=1.0e-7;
crit=0.0004;
NCini=5; %initially check five points
while (critv>diffv)
    vrfunc_old=vrfunc;
          for ia=1:NA
            if ia==1
                cmin=climit;
            else
            cmin=crfunc(ia-1);
            end
            cmax=(eJr+(1+r)*Avector(ia)-tax_T3(eJr,r*Avector(ia)));
            cg=linspace(cmin,cmax,NCini);
            afuture=eJr+(1+r)*Avector(ia)-tax_T3(eJr,r*Avector(ia))-cg;
            vtemp=uc(cg)+b*s*interp1(Avector,vrfunc,min(amax,afuture),'spline');
            [vtemp_v, ctemp_g]=max(vtemp);
            cmax_v=cg(ctemp_g);
            vmax_v=vtemp_v;
            
            diff=(cmax-cmin)/4;
            while diff>crit

            diff=diff/2;
          
            if ((cmax_v>cmin))
                cg=cmax_v-diff;
                afuture=eJr+(1+r)*Avector(ia)-tax_T3(eJr,r*Avector(ia))-cg;
                vtemp=uc(cg)+b*s*interp1(Avector,vrfunc,min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                    continue
                end %if
            end %while internal
            if ((cmax_v<cmax))
                cg=cmax_v+diff;
                afuture=eJr+(1+r)*Avector(ia)-tax_T3(eJr,r*Avector(ia))-cg;
                vtemp=uc(cg)+b*s*interp1(Avector,vrfunc,min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                end %if
            end %while internal
            end %while
            vrfunc(ia)=vmax_v;
            crfunc(ia)=cmax_v;
     
           end %ia
            
critv=max(abs(vrfunc_old-vrfunc));
end

vfunc_temp=zeros(NA,1);
cfunc_temp=zeros(NA,1);
vfunc{JR}=vrfunc;
cfunc{JR}=crfunc;
for j=JR-1:-1:1 %retirees' life
    vfuture=vfunc{j+1};
        for ia=1:NA
            if ia==1
                cmin=climit;
            else
            cmin=cfunc_temp(ia-1);
            end
            cmax=(e(j)+(1+r)*Avector(ia)-tax_T3(e(j),r*Avector(ia)));
            cg=linspace(cmin,cmax,NCini);
            afuture=e(j)+(1+r)*Avector(ia)-tax_T3(e(j),r*Avector(ia))-cg;
            vtemp=uc(cg)+b*interp1(Avector,vfuture,min(amax,afuture),'spline');
            [vtemp_v, ctemp_g]=max(vtemp);
            cmax_v=cg(ctemp_g);
            vmax_v=vtemp_v;
            
            diff=(cmax-cmin)/4;
            while diff>crit

            diff=diff/2;
          
            if ((cmax_v>cmin))
                cg=cmax_v-diff;
                afuture=e(j)+(1+r)*Avector(ia)-tax_T3(e(j),r*Avector(ia))-cg;
                vtemp=uc(cg)+b*interp1(Avector,vfuture,min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                    continue
                end %if
            end %while internal
            if ((cmax_v<cmax))
                cg=cmax_v+diff;
                afuture=e(j)+(1+r)*Avector(ia)-tax_T3(e(j),r*Avector(ia))-cg;
                vtemp=uc(cg)+b*interp1(Avector,vfuture,min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                end %if
            end %while internal
            end %while
            vfunc_temp(ia)=vmax_v;
            cfunc_temp(ia)=cmax_v;
            %vfunc{j}(ia)=vmax_v;
            %cfunc{j}(ia)=cmax_v;
     
        end %ia
        vfunc{j}=vfunc_temp;
        cfunc{j}=cfunc_temp;
end %j


%=====================
%simulate optimal decisions
%=====================
JN=100; %maximum periods
vopt=zeros(JN,1); %optimal value
copt=zeros(JN,1); %consumption
aopt=zeros(JN+1,1); %asset
taxopt=zeros(JN,1); %tax payments
taxableopt=zeros(JN,1); %taxable income
for j=1:JR-1
            vfuture=vfunc{j+1};
            cmin=climit;
            cmax=(e(j)+(1+r)*aopt(j)-tax_T3(e(j),r*aopt(j)));
            cg=linspace(cmin,cmax,NCini);
            afuture=e(j)+(1+r)*aopt(j)-tax_T3(e(j),r*aopt(j))-cg;
            vtemp=uc(cg)+b*interp1(Avector,vfuture, min(amax,afuture),'spline');
            [vtemp_v, ctemp_g]=max(vtemp);
            cmax_v=cg(ctemp_g);
            vmax_v=vtemp_v;
            
            diff=(cmax-cmin)/4;
            while diff>crit

            diff=diff/2;
          
            if ((cmax_v>cmin))
                cg=cmax_v-diff;
                afuture=e(j)+(1+r)*aopt(j)-tax_T3(e(j),r*aopt(j))-cg;
                vtemp=uc(cg)+b*interp1(Avector,vfuture, min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                    continue
                end %if
            end %while internal
            if ((cmax_v<cmax))
                cg=cmax_v+diff;
                afuture=e(j)+(1+r)*aopt(j)-tax_T3(e(j),r*aopt(j))-cg;
                vtemp=uc(cg)+b*interp1(Avector,vfuture, min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                end %if
            end %while internal
            end %while
            copt(j)=min(round(cmax_v,3),floor((e(j)+(1+r)*aopt(j)-tax_T3(e(j),r*aopt(j)))*1000)/1000);          

    aopt(j+1)=e(j)+(1+r)*aopt(j)-tax_T3(e(j),r*aopt(j))-copt(j);   
    taxopt(j)=tax_T3(e(j),r*aopt(j));
    taxableopt(j)=e(j)+r*aopt(j);
end

for j=JR:JN 
            vfuture=vrfunc;
            cmin=climit;
            cmax=(eJr+(1+r)*aopt(j)-tax_T3(eJr,r*aopt(j)));
            cg=linspace(cmin,cmax,NCini);
            afuture=eJr+(1+r)*aopt(j)-tax_T3(eJr,r*aopt(j))-cg;
            vtemp=uc(cg)+b*s*interp1(Avector,vfuture, min(amax,afuture),'spline');
            [vtemp_v, ctemp_g]=max(vtemp);
            cmax_v=cg(ctemp_g);
            vmax_v=vtemp_v;
            
            diff=(cmax-cmin)/4;
            while diff>crit

            diff=diff/2;
          
            if ((cmax_v>cmin))
                cg=cmax_v-diff;
                afuture=eJr+(1+r)*aopt(j)-tax_T3(eJr,r*aopt(j))-cg;
                vtemp=uc(cg)+b*s*interp1(Avector,vfuture, min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                    continue
                end %if
            end %while internal
            if ((cmax_v<cmax))
                cg=cmax_v+diff;
                afuture=eJr+(1+r)*aopt(j)-tax_T3(eJr,r*aopt(j))-cg;
                vtemp=uc(cg)+b*s*interp1(Avector,vfuture, min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                end %if
            end %while internal
            end %while
            copt(j)=min(round(cmax_v,3),floor((eJr+(1+r)*aopt(j)-tax_T3(eJr,r*aopt(j)))*1000)/1000);          

    aopt(j+1)=eJr+(1+r)*aopt(j)-tax_T3(eJr,r*aopt(j))-copt(j);   
    taxopt(j)=tax_T3(eJr,r*aopt(j));
    taxableopt(j)=eJr+r*aopt(j);
end


%cumulative utility if adopting the optimal strategy
paymentopt=zeros(JN,1);
j=1;
paymentopt(j)=b^(j-1)*uc(copt(j));
for j=2:JN
paymentopt(j)=paymentopt(j-1)+b^(j-1)*uc(copt(j));
end




%hand to mouth consumers
paymenthandtomouth=zeros(JN,1);
j=1;
paymenthandtomouth(j)=b^(j-1)*uc(e(j)-tax_T3(e(j),0));
for j=2:JR-1
paymenthandtomouth(j)=paymenthandtomouth(j-1)+b^(j-1)*uc(e(j)-tax_T3(e(j),0));
end
for j=JR:JN
paymenthandtomouth(j)=paymenthandtomouth(j-1)+b^(j-1)*uc(eJr-tax_T3(eJr,0));
end

chandtomouth=zeros(JN,1);
j=1;
chandtomouth(j)=e(j)-tax_T3(e(j),0);
for j=2:JR-1
chandtomouth(j)=(e(j)-tax_T3(e(j),0));
end
for j=JR:JN
chandtomouth(j)=(eJr-tax_T3(eJr,0));
end

prob_d=zeros(JN+1,1); %the probability of decease in period j
for j=JR+1:JN
    prob_d(j)=prob_d(j-1)+(1-s)*(1-prob_d(j-1));
end
prob_d(JN+1)=1; %maximum age of JN




copt_T3=copt;
paymentopt_T3=paymentopt;
taxableopt_T3=taxableopt;
taxopt_T3=taxopt;
sopt_T3=zeros(JN,1);
aopt_T3=aopt;


for j=1:JR-1
sopt_T3(j)=e(j)-taxopt_T3(j)-copt_T3(j);
end

for j=JR:JN
sopt_T3(j)=eJr-taxopt_T3(j)-copt_T3(j);    
end


fprintf('T3: Selected aggregate statistics')
[sum((prob_d(2:JN+1)-prob_d(1:JN)).*paymentopt_T3); paymentopt_T3(15); sum((prob_d(2:JN+1)-prob_d(1:JN)).*paymentopt_T3)-paymentopt_T3(15);...
    0; aopt_T3(JR);0;  mean(copt_T3(1:JR-1));  mean(copt_T3(16:25)); copt_T3(16); sum((1-prob_d(1:JN)).*taxopt_T3); sum(sopt_T3(1:JR-1))]

save T3opt.mat copt_T3 paymentopt_T3 sopt_T3  aopt_T3 taxableopt_T3 taxopt_T3 z  prob_d


fprintf('Treatment 3\n')
fprintf('1.consumption   2.cumulative earnings   3.saving    4.asset balance  5.ordinary income   6.taxes  7.survival prob.')
T3=[copt_T3(1:JN) paymentopt_T3(1:JN) sopt_T3(1:JN)  aopt_T3(1:JN) taxableopt_T3(1:JN) taxopt_T3(1:JN) (1-prob_d(1:JN))]






%=====================
%Appendix Figure C1 
%=====================

Jgraph=35; %only show 30 periods for graphs
close all
fName = 'Times';
fSize = 8;
 
%   Paper size and position for single figures
pSize1     = [4.2,3.2];
pPosition1 = [0.1,0.1,4.0,3.0];
 
%   Paper size and position for double figures
pSize2     = [3.2,3.0];
pPosition2 = [0,0,3.2,3.0];

figure

axes('YGrid','off','XGrid','off',...
     'FontName',fName,'FontSize',fSize);
hold on
box on
scatter(eJr-tax(eJr),uc(eJr-tax(eJr)),10,[0 0 0])
plot(linspace(47.4,200,100),uc(linspace(47.4,200,100)),'LineWidth',1,'Color',[0 0 0])
yl=ylabel('Utility/ Payoff');
xl=xlabel('Consumption (k)');
lg=legend('Retirement consumption with no savings');
set(lg,'Location','SouthEast','FontName',fName,'FontSize',fSize)
set(xl,'FontName',fName,'FontSize',fSize);
set(yl,'FontName',fName,'FontSize',fSize);
set(gca,'FontName',fName,'FontSize',fSize);
set(gcf,'paperunits','inches')
set(gcf,'papersize',pSize1)
set(gcf,'PaperPosition',pPosition1)
print('-dpdf','-r600','graph\utility.pdf')


figure
axes('YGrid','off','XGrid','off',...
     'FontName',fName,'FontSize',fSize);
hold on
box on
plot(linspace(0,200,100),tax(linspace(0,200,100)),'LineWidth',1,'Color',[0 0 0])
yl=ylabel('Tax (k)');
xl=xlabel('Taxable income (k)');
set(xl,'FontName',fName,'FontSize',fSize);
set(yl,'FontName',fName,'FontSize',fSize);
set(gca,'FontName',fName,'FontSize',fSize);
set(gcf,'paperunits','inches')
set(gcf,'papersize',pSize1)
set(gcf,'PaperPosition',pPosition1)
print('-dpdf','-r600','graph\taxfunction.pdf')


figure
axes('YGrid','off','XGrid','off',...
     'FontName',fName,'FontSize',fSize);
hold on
box on
plot(1:Jgraph,[e; repmat(eJr,Jgraph-JR+1,1)],'LineWidth',1,'Color',[0 0 0])
line([JR JR],[0 max(e)],'LineWidth',1,'Color','Blue')
yl=ylabel('Endowment (k)');
xl=xlabel('Period');
set(xl,'FontName',fName,'FontSize',fSize);
set(yl,'FontName',fName,'FontSize',fSize);
set(gca,'FontName',fName,'FontSize',fSize);
set(gcf,'paperunits','inches')
set(gcf,'papersize',pSize1)
set(gcf,'PaperPosition',pPosition1)
print('-dpdf','-r600','graph\endowment.pdf')


figure
axes('YGrid','off','XGrid','off',...
     'FontName',fName,'FontSize',fSize);
hold on
box on
plot(1:Jgraph,s.^(max(0,(1:Jgraph)-16)),'LineWidth',1,'Color',[0 0 0])
line([JR JR],[0 1],'LineWidth',1,'Color','Blue')
yl=ylabel('Unconditional survival rate');
xl=xlabel('Period');
set(xl,'FontName',fName,'FontSize',fSize);
set(yl,'FontName',fName,'FontSize',fSize);
set(gca,'FontName',fName,'FontSize',fSize);
set(gcf,'paperunits','inches')
set(gcf,'papersize',pSize1)
set(gcf,'PaperPosition',pPosition1)
print('-dpdf','-r600','graph\s.pdf')



